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A stochastic leap-frog algorithm for the numerical integration of Brownian motion stochastic differ- 
ential equations with multiplicative noise is proposed and tested. The algorithm has a second-order 
convergence of moments in a finite time interval and requires the sampling of only one uniformly dis- 
tributed random variable per time step. The noise may be white or colored. We apply the algorithm 
to a study of the approach towards equilibrium of an oscillator coupled nonlinearly to a heat bath 
and investigate the effect of the multiplicative noise (arising from the nonlinear coupling) on the 
relaxation time. This allows us to test the regime of validity of the energy-envelope approximation 
method. 
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I. INTRODUCTION 

Stochastic differential equations with multiplicative 
noise have not only found many applications in physics 
but also have interesting mathematical properties. Con- 
sequently they have attracted substantial attention over 
the years ]j|-|lT|]. The key point lies in the fundamen- 
tal difference between additive and multiplicative noises: 
Additive noise does not couple directly to the system 
variables and disappears from the noise-averaged form of 
the dynamical equations. However, in the case of multi- 
plicative noise, the system variables do couple directly to 
the noise (alternatively, we may say that the noise am- 
plitude depends on the system variables). This fact can 
lead to dramatic changes of system behavior that cannot 
occur in the presence of additive noise alone. Two classic 
illustrations are the Kubo oscillator |12| and the existence 



of long-time tails in transport theory |13| . In this paper 
we will investigate another example, that of an oscillator 
nonlinearly coupled to a heat bath, in which the effects 
of multiplicative noise can significantly alter the qualita- 
tive nature, as well as the rate S, of the equilibration 
process (relative to that of an oscillator subjected only 
to additive noise). 

The dynamical behavior of systems subjected to noise 
can be studied in two different ways: we may either solve 
stochastic differential equations and average over realiza- 
tions to obtain statistical information, or we may directly 
solve the Fokker-Planck equation which describes the 
evolution of the corresponding probability distribution 
function. Both approaches have their share of advantages 
and disadvantages. Fokker-Planck equations are partial 
differential equations and their mathematical properties 
are still not fully understood. Moreover, they are very 
expensive to solve numerically even for dynamical sys- 
tems possessing only a very modest number of degrees 
of freedom. Truncation schemes or closures (such as cu- 
mulant truncations) have had some success in extracting 



the behavior of low-order moments, but the systematics 
of these approximations remains to be elucidated. Com- 
pared to the Fokker-Planck equation, stochastic differ- 
ential equations are not difficult to solve, and with the 
advent of modern supercomputers, it is possible to run 
very large numbers of realizations in order to compute 
low-order moments accurately. (We may mention that in 
applications to field theories it is essentially impossible to 
solve the corresponding Fokker-Planck equation since the 
probability distribution is now a functional.) However, 
the extraction of the probability distribution function it- 
self is very difficult due to the sampling noise inherent in 
a particle representation of a smooth distribution. 

Numerical algorithms to solve stochastic differential 
equations have been discussed extensively in the litera- 
ture 14 - 19 . The simplest, fastest, and still widely-used, 
is Euler's method which yields first-order convergence of 
moments for a finite time interval. Depending on the 
control over statistical errors arising from the necessarily 
finite number of realizations, in the extraction of statis- 
tical information it may or may not pay to use a higher 
order algorithm especially if it is computationally expen- 
sive. Because of this fact, it is rare to find high-order 
schemes being put to practical use for the solution of 
stochastic differential equations, and second-order con- 
vergence is usually considered a good compromise be- 
tween efficiency and accuracy. A popular algorithm with 
second-order convergence of moments for additive noise 
but with only first-order convergence of moments for mul- 
tiplicative noise is Heun's algorithm (also called stochas- 
tic RK2 by some authors) [ pd|JT7|j20[| . A stochastic leap- 
frog algorithm which has the same order convergence of 
moments as Heun's method was suggested in Ref. 21 1 to 
study particle motion in a stochastic potential without 
damping. Several other algorithms for particle motion in 
a quasi-conservative stochastic system were proposed in 
Ref. JT|] and in the book by Allen and Tildesley ||. At 
every time step, these methods all require sampling two 
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Gaussian random variables which adds to the computa- 
tional cost. A modified algorithm suggested in Ref. jllj 
requires only one Gaussian random variable but applies 
only to white Gaussian noise. In the following sections, 
we present a new stochastic leap-frog algorithm for mul- 
tiplicative Gaussian white noise and Ornstein-Uhlenbeck 
colored noise which not only has second-order conver- 
gence of moments but also requires the sampling of only 
one random uniform variable per time step. 

The organization of this paper is as follows: General 
numerical integration of a system of stochastic differ- 
ential equations with Gaussian white noise is discussed 
in Section II. The stochastic leap-frog algorithms for 
Brownian motion with Gaussian white noise and colored 
Ornstein-Uhlenbeck noise are given in Section III. Nu- 
merical tests of these algorithms using a one-dimensional 
harmonic oscillator are presented in Section IV. A phys- 
ical application of the algorithm to the multiplicative- 
noise Brownian oscillator is given in Section V. Section VI 
contains the final conclusions and and a short discussion. 



II. NUMERICAL INTEGRATION OF 
STOCHASTIC DIFFERENTIAL EQUATIONS 

A general system of continuous- time stochastic differ- 
ential equations (Langevin equations) can be written as 
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where i — 1, • ■ • , n and £j(t) is a Gaussian white noise 
with 
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and the symbol (■ ■ •) represents an average over realiza- 
tions of the inscribed variable (ensemble average). The 
noise is said to be additive when <jy is not a function 
of the Xi, otherwise it is said to be multiplicative. In 
the case of multiplicative noises, a mathematical subtlety 
arises in interpreting stochastic integrals, the so-called 
Ito-Stratonovich ambiguity p3[ . It should be stressed 
that this is a point of mathematics and not of physics. 
Once it is clear how a particular Langevin equation has 
been derived and what it is supposed to represent, it 
should either be free of this ambiguity (as in the case of 
the example we study later) or it should be clear that 
there must exist two different stochastic equations, one 
written in the Ito form, the other in Stratonovich, both 
representing the same physical process and hence yielding 
identical answers for the variables of interest. (Another 
way to state this is that there should be only one unique 
Fokker-Planck equation.) It is important to note that the 
vast majority of numerical update schemes for Langevin 
equations use the Ito form of the equation. 

The integral representation of the set of equations (|l|) 

is 



Xi(t) = Xi(0) + / dsF i (x 1 (s), - ■ ■ ,x n (s)) 
Jo 

+ dsaij(xi(s),- ■ ■ ,x n (s))£j(s) (4) 
Jo 

where Xi(0) is a given sharp initial condition at t = 0. 
The infinitesimal update form of this equation may be 
derived by replacing t with an infinitesimal time step ft: 
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Since Fi and try are smooth functions of the Xi, they may 
be expanded about their values at t — 0, in which case 
we can write the exact solution for x% (ft) as 



Xi {h) = Di(h) + Si(h) 



(6) 



where Di(h) and Si (ft) denote the deterministic and 
stochastic contributions respectively. The deterministic 
contribution Di{h) is 

A(ft) = a!i(0) + hFi + h 2 F hk F k + 0(h 3 ) (7) 

where i 7 ^ = dFi/dx k , the summation convention for the 
repeated indices having being employed. The stochastic 
contribution Si{h) is 

S,(ft) = a i0 Wj(h) + <T ljtk a k iCij(h) + F lik a kl Zi(h) 

+<r ij , k F k (hW j (h) - Zj (h)) + l ri , kl a km a ln H mn] {h) 

+ -^F i ^ k ia ks ai t G st (h) + -F k a i3 , k i(7i m K m j(h) 
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The quantities W%, Cy, H ljkl Z i} G ljl K tj , and I i]kl 
are random variables which can be written as stochastic 
integrals over the Gaussian white noise £(i): 



Wi(h) = \ dt&(t) ~ Oih 1 ' 2 ) 



(9) 
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H ijk (h)= / dtWiWWjft&it) ~ 0(h 3 ' 2 ) (11) 
Jo 



Cy(ft)= / dtWi(t)^(t) ~ 0(h) 
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Zi{h) = / dtWi(t) ~ 0(/i 3 / 2 ) (12) 
Jo 
r h 

Gij(h)= dtWi(t)Wj(t) ~ 0(h 2 ) (13) 
Jo 

rh 

Kij(h)= / tdtWiW&it) ~ 0(h 2 ) (14) 
Jo 

rh 

I ijk i{h) = dtW^W^Wkitfoit) ~ 0(h 2 ) (15) 



Ito integration has been employed in the derivation of 
the above equations. 

The nth moment of the xi is 

(xi(h) n ) = ((Diih) + Si(h)) n ) 

= D, l (h) n + nD l {h) n - 1 (S t (h)) 

+ClD i {h) n -*\{S i {h)f)+--- (16) 



where 



CI 



n 



and 



n J i\(n - i)\ 
1 



(17) 



(s l (h)) = l F: k i° ks v lSh2 +°( h3 ) 

(Si(h)S:j(h)) = a a a jl h + ^Ja kl (j^a pl h 2 

+ \<? U F j k a kl h 2 + ^a jl F* k a kl h 2 



1 



1 



+ -a il a j lF k h 2 + -a jl a%F k h 2 
+^a ip a%a km a lm h 2 

+\v jp v%i<r km ° lm h 2 + 0(h 3 ) (19) 

{S i (h)S j (h)S k (h)) = 0(h 3 ) (20) 
(S^h) 4 ) = 3(a u ) 4 + 0(h 3 ) (21) 
((S i (h)f)=0(h s ) (22) 

Suppose that the results from a numerical algorithm were 
written as 



Xi(h) = Di(h) + Si(h) 



(23) 



where the Xi are approximations to the exact solutions 
Xi . The nth moment of Xi is 

(xi(h) n ) = {(Di(h) + Si(h)) n ) 

= D t (h) n + nD i {h) n ~ 1 (S t (h)) 

+C 2 l D l (hr~ 2 ((S l (h)) 2 } + --- (24) 

Comparing Eqns. ( |l6| ) and (p4|), we see that if Di(h) and 
Di(h), and and coincide up to ft, 2 , we will 

have 



Xi (h) - Xi{h) = 0(h 3 ) 
and for a finite time interval 

(x i (tr)-{x i (t)) n ) = o(h 2 ) 



(25) 
(26) 



III. STOCHASTIC LEAP-FROG ALGORITHM 
FOR BROWNIAN MOTION 

The approach to modeling Brownian motion that we 
consider here is that of a particle coupled to the environ- 
ment through its position variable Q. When this is the 
case, noise terms enter only in the dynamical equations 
for the particle momenta. In the case of three dimensions, 
the dynamical equations take the general form: 

ii = F 1 (X 1 ,X2,X 3 ,X 4: ,X 5 ,X 6 ) + <7 11 (x 2l X i ,XQ)^ 1 (t) 

x 2 = F 2 (xt) 

*3 = F 3 (xi,x 2 ,x 3 ,x 4:1 x 5 ,x 6 ) + cr 33 (x 2l x 4: ,x 6 )^ 3 (t) 
±4 = F 4 (x 3 ) 

±5 = F 5 (xi,x 2 ,x 3 ,x 4:1 x 5 ,x 6 ) + cr 55 (x 2l x 4: ,x 6 )^ 5 {t) 
xe = F e (x 5 ) (27) 

The convention used here is that the odd indices corre- 
spond to momenta, and the even indices to the spatial 
coordinate. In the dynamical equations for the momenta, 
the first term on the right hand side is a systematic drift 
term which includes the effects due to external forces and 
damping. The second term is stochastic in nature and 
describes a noise force which, in general, is a function of 
position. The noise £(i) is first assumed to be Gaussian 
and white as defined by Eqns. 
leap-frog algorithm for the Eqns 




). The stochastic 
7|) is written as 



Xi(h) = Di(h)+Si(h) 



(28) 



The deterministic contribution Di(h) can be obtained us- 
ing the deterministic leap-frog algorithm. The stochastic 
contribution Si(h) can be obtained by applying Eq. (JsJ) on 
Eq. (|27j ). The stochastic integration defined by Eqs. (||) 
to (|15|) can be approximated so that the moment rela- 
tionships defined by Eqs. ( |l8| ) to (|2^) are satisfied. After 
some calculation, the deterministic contribution Di(h) 
and the stochastic contribution Si{h) of the above recur- 
sion formula for one-step integration are found to be 

Di(h) = Xi(0) + hF^xl, X2, xl, x* 4 , xl, xl); 

{i = 1,3,5} 
Di(h) = x* 

-\-—hFi \xi-\ + /i-Fj_i (x^ , x 2 , -^3? ®5) ^e)J ! 
{i = 2,4,6} 
Si{h) = a u VhWi(h) + ^F i<k u kk h^ 2 Wi{h) 

+\v ii>j F j h?' 2 W i {h) 

+ ^F hkl a kk cr ll h 2 W l (h)W l (h); 
{i = 1,3,5; j = 2,4,6; k,l= 1,3,5} 
Si(h) = -^F^h^W^h) 
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+ -F hJJ a 2 JJ h 2 W J (h)W J (h) 

{* = 2,4,6; j = 1,3,5} 

sci(0) + -hF^xi, x 2 ,x 3 ,X4,x 5 ,x e ) 

= 1,2,3,4,5,6} 



(29) 



where Wi(ft) is a series of random numbers with the mo- 
ments 

(W t (h)) = ((Wi(h)) 3 ) = (m(h)f) = (30) 

m(h)) 2 ) = i, ((^w) 4 ) = 3 (si) 

This can not only be achieved by choosing true Gaus- 
sian random numbers, but also by using the sequence of 
random numbers following: 



-V3, R < 1/6 
Wi(h) = { 0, 1/6 < R < 5/6 
V3, 5/6 < i? 



(32) 



where R is a uniformly distributed random number on 
the interval (0,1). This trick significantly reduces the 
computational cost in generating random numbers. 

Next we consider the case that the noise in Eqs. ( pjj ) 
is a colored Ornstein-Uhlenbeck process which obeys 



Si(h) = -^=a l i(x2,X4,x 6 )k i h 3/ ' 2 Wi(h)- 1 
V3 

0=1,3,5} 
Si(h) = 0; 

{i = 2,4,6} 

4 = kiVhWi(h) - h 2 h 3 ^Wi(h); 



where 



= 1,3,5} (36) 



x* = Xi(Q) + -h [F i (x 1 ,x 2 ,x 3 ,x 4 ,x 5 ,x 6 ) 

+o-ii(x2,X4,x 6 )£i} ; 
{i = 1,3,5} 

x* = Xi(0) + -hFi(xi, X2, x 3 , x 4 , x 5 , xe); 
{i = 2,4,6} 

e? = €i(0)exp(-ifci/i)i 

= 1,3,5} (37) 



<6(*)> = 



(6lW)> 



exp(-fci|t-t'|) 



(33) 
(34) 



where the correlation factor k% is the reciprocal of the 
correlation time. In the limit of fej — * oo, the Ornstein- 
Uhlenbeck process reduces to Gaussian white noise. The 
above process can be generated by using a white Gaussian 
noise from a stochastic differential equation 



(35) 



where d(t) is a standard Gaussian white noise. The ini- 
tial value £j (0) is chosen to be a Gaussian random number 
with (6(0)) = and (&(0) 2 ) = 

For the stochastic process with colored noise, the leap- 
frog algorithm for Eqns. ( P7|) is of the same form as that 
for white noise (Cf. Eqn. (|29|)), but with 

Di(h) = Xi(0) + hFi{x~l,x^,xl,x%,xl,x%) 
-\-hau (x 2 , X4 , x e ; 
= 1,3,5} 
x* 



Di(h) 



D 6i (h) 



+~hFi [xi-i + hFi-^xl, x 2 * , x* 3 , x A * , x* 5 , x 6 *) 

+hcri-u-i(xZ,xl,xl)($_ 1 ] ; 
= 2,4,6} 
£j(0)exp(-fcj/i); 
0=1,3,5} 



FIG. 1. Zero damping convergence test. Top: (x 2 (t)) at 
t = 6 as a function of step size with white Gaussian noise. 
Bottom: (x 2 (t)) at t — 6 as a function of step size with colored 
Ornstein-Uhlenbeck noise. Solid lines represent quadratic fits 
to the data points (diamonds). 



IV. NUMERICAL TESTS 

The above algorithms were tested on a one-dimensional 
stochastic harmonic oscillator with a simple form of the 
multiplicative noise. The equations of motion were 
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p = F 1 {p,x) + a(x)at) 
x = p 



(38) 



where F\(p,x) — —~/p— rj 2 x and <r(x) — —ax. 

As a first test, we computed (a; 2 ) as a function of time 
step size. To begin, we took the case of zero damping 
constant (7 = 0), where (a; 2 ) can be determined analyt- 
ically. The top curve in Fig. [l] shows (a; 2 ) at t = 6.0 as 
a function of time step size with white Gaussian noise. 
Here, the parameters 7/ and a are set to 1.0 and 0.1. 
The ensemble averages were taken over 10 6 independent 
simulations. The analytically determined value of (x 2 ) 
at t = 6.0 is 2.095222 (The derivation of the analytical 
results is given in the Appendix). The quadratic con- 
vergence of the stochastic leap-frog algorithm is clearly 
seen in the numerical results. We then considered the 
case of colored Ornstein-Uhlenbeck noise as a function of 
time step size using the same parameters as in the white 
Gaussian noise case and with the correlation parameter 
k = 0.16. The result is shown as the bottom curve in 
Fig. [l] and the quadratic convergence is again apparent. 



gorithm, we computed (x 2 ) as a function of t using 
100, 000 numerical realizations for a particle starting 
from (0.0, 1.5) in the {x,p) phase space. The results along 
with the analytical solution and a numerical solution us- 
ing Heun's algorithm are given in Fig. Parameters 
used were h = 0.1, 77 = 1.0, and a — 0.1. The ad- 
vantage in accuracy of the stochastic leap-frog algorithm 
over Heun's algorithm is clearly displayed, both in terms 
of error amplitude and lack of a systematic drift. 

We note that while in general Heun's algorithm is only 
linear for multiplicative noise applications, for the partic- 
ular problem at hand it turns out to be quadratic. This 
is due to a coincidence: the stochastic term of x does 
not contain W(h) but does posses a higher order term 
hW(h). However, this higher order term has a larger 
coefficient compared with our stochastic leap-frog algo- 
rithm, and this accounts for the larger errors observed in 
Fig. 3. 
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FIG. 3. Comparing stochastic leap-frog and the Heun al- 
gorithm: (x 2 (t)} as a function of t. Errors are given relative 
to the exact solution. 



FIG. 2. Finite damping (7 = 0.1) convergence test. Top: 
(x 2 {t)} a,tt = 12 as a function of step size with white Gaussian 
noise. Bottom: (x 2 (t)} at t = 12 as a function of step size 
with colored Ornstein-Uhlenbeck noise. Solid lines represent 
quadratic fits to the data points (diamonds). 

We verified that the quadratic convergence is present 
for nonzero damping (7 = 0.1). At t — 12.0, and with 
all other parameters as above, the convergence of (x 2 ) 
as a function of time step is shown by the top and bot- 
tom curves in Fig. ^ (white Gaussian noise and colored 
Ornstein-Uhlenbeck noise, respectively). 

As a comparison against the conventional Heun's al- 



V. A PHYSICAL APPLICATION: THE 
MECHANICAL OSCILLATOR 

In this section, we apply our algorithm to studying the 
approach to thermal equilibrium of an oscillator coupled 
nonlinearly to a heat bath modeled by a set of nonintcr- 
acting harmonic oscillators jjj]. The nonlinear coupling 
leads to the introduction of multiplicative noise into the 
system dynamics. Lindenberg and Seshadri have pointed 
out that, at weak coupling, multiplicative noise may sig- 
nificantly enhance the equilibration rate relative to the 
rate for weak linear coupling (additive noise) |2j . We will 
choose the same form of the coordinate couplings as in 
Ref. 0, in which case the additive noise equations are 



p = -u%x - X p + y/2D £ (t) 
x = p 



(39) 
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and for the system with multiplicative noise only: 



p = —UIqX 

x = p 



-~-\ 2 x 2 p-^2D a x&(t) 



(40) 



where the diffusion coefficients Di = X{kT, i = 0, 2, A$ is 
the coupling constant, k is Boltzmann's constant, T is the 
heat bath temperature, and luq is the oscillator angular 
frequency without damping. The approach to thermal 
equilibrium is guaranteed for both sorts of noises by the 
fluctuation-dissipation relation 



(41) 



written here for the general case when both noises are 
simultaneously present. While in all cases, it is clear 
that the final distribution is identical and has to be the 
thermal distribution, the precise nature of the approach 
to equilibrium can certainly be different. We wish to 
explore this issue in more detail. An important point to 
keep in mind is that in this particular system of equations 
there is no noise-induced drift in the Fokker-Planck equa- 
tion obtained from the Stratonovich form of the Langevin 
equation, i.e., there is no Ito-Stratonovich ambiguity. 

It is a simple matter to solve the Langevin equations 
given above applying the algorithm from Eqs. (p9|) . As 
our primary diagnostic, we computed the noise-averaged 
energy (E(t)) of the oscillator as a function of time t, 
where 



E(t) = \p 



2,1 2 2 

2^ ■ 



(42) 



In the weak coupling limit and employing orbit-averaging 
(valid presumably when the dynamical time scale is much 
smaller than the relaxation time scale) , one finds [0 



(E(t)) =kT- (kT-E )e 



-Ant 



(43) 



in the case of additive noise (a result which can also be 
directly obtained as a limiting case from the known form 
of the exact solution given, e.g., in Ref. J24|). The corre- 
sponding form of the approximate solution in the case of 
multiplicative noise is 



(E(t)) 



E kT 



E Q + (kT - E Q ) cxp(-A 2 m/cj2) ' 



(44) 



While in the case of additive noise, the exponential na- 
ture of the relaxation is already clear from the form of the 
exact solution (cf. Ref. 0]), the situation in the case of 
multiplicative noise is not obviously apparent as no exact 
solution is known to exist. The prediction of a relaxation 
process controlled by a single exponential as found in ([44| ) 
is a consequence of the assumption (x 2 (t)) ~ kT/w 2 at 
"late" times, this implying a constant damping coefficient 
in the Langevin equation (|40|). 

The timescale separations necessary for the energy- 
envelope method to be applicable are encoded in the fol- 
lowing inequalities gj: 



kTX 2 

— 5- < l; 



additive noise 



multiplicative noise 



(45) 
(46) 



As a first check, we performed simulations with ujq = 1.0, 
Ao = A2 = 0.01, and kT — 4.5, in which case both 
the above conditions are satisfied. Moreover, with these 
choices of parameter values, and within the energy en- 
velope approximation, the relaxation time predicted for 
multiplicative noise is substantially smaller than for the 
case of additive noise. At the same time we also ran a 
simulation at kT = 200 to see how the energy envelope 
approximation for multiplicative noise breaks down at 
high temperatures. 




FIG. 4. Temporal evolution of the scaled average energy 
{E(t))/kT with additive noise and multiplicative noise. The 
dashed lines I and II are the predictions from Eqn. (^ij) for 
kT = 200 and kT = 4.5 respectively. The dashed line III is 
the theoretical prediction for additive noise with kT = 4.5. 
As predicted, the relaxation proceeds much faster with multi- 
plicative noise: The solid lines are numerical results for mul- 
tiplicative noise at kT = 200 and kT — 4.5. It is clear that at 
higher temperatures, the theory grossly underestimates the 
relaxation time. 

In Fig. ^, we display the time evolution of the aver- 
age energy (scaled by kT for convenience) with additive 
and multiplicative noise both from the simulations and 
the approximate analytical calculations. In the case of 
weak coupling to the environment (small Ao, A2), the 
rate at which the average energy approaches equilibrium 
is significantly greater for the case of multiplicative noise 
relative to the case of additive noise more or less as ex- 
pected. In addition, the analytic approximation result- 
ing from the application of the energy-envelope method 
( |44| ) is seen to be in reasonable agreement with the nu- 
merical simulations for kT = 4.5. The slightly higher 
cqulibration rate from the analytical calculation is due to 
the truncation in the energy envelope equation using the 
(E 2 (t)} re 2(E(t)} 2 relation which yields an upper bound 



G 



on the rate of equilibration of the average energy [Q. 
Note that in the case of high temperature (kT = 200) 
the relaxaton time computed from the energy envelope 
method is much smaller than the numerical result, con- 
sistent with the violation of the condition (ff6|). 

While the results shown in Fig. [| do show that the 
energy envelope approximation is qualitatively correct 
within its putative domain of validity, it is clear that 
the actual relaxation process is not of the precise form 
(fl4|). In Fig. ||we illustrate this point by plotting 



EpjkT - (E(t))) 
(E(t))(kT - E ) 



= exp(— \2kTt / ujq) 



(47) 



[equivalent to (44)] against time on a log scale: the re- 
laxation is clearly nonexponential. The reason for the 
failure of the approximation is that despite the fact that 
equipartition of energy does take place on a relatively 
short time scale, it is not true that (x 2 (t)) can be treated 
as a constant even at relatively late times. 



oscillator-heat-bath system in order to investigate the ef- 
fect of multiplicative noise on the nature of the relaxation 
process. 
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FIG. 5. The LHS of @ as a function of time (straight line) 
compared with numerical results for kT = 4.5. Also shown 
is a numerical result for the case of additive noise which is in 
excellent agreement with the predicted exponential relaxation 
with the relaxation timescale = 1/Aq. 

VI. CONCLUSIONS 

We have presented a stochastic leap-frog algorithm 
for single particle Brownian motion with multiplicative 
noise. This method has the advantages of retaining the 
symplcctic property in the deterministic limit, ease of im- 
plementation, and second-order convergence of moments 
for multiplicative noise. Sampling a uniform distribution 
instead of a Gaussian distribution helps to significantly 
reduce the computational cost. A comparison with the 
conventional Heun's algorithm highlights the gain in acu- 
racy due to the new method. Finally, we have applied the 
stochastic leap-frog algorithm to a nonlinearly coupled 
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APPENDIX A: 

The analytic solution of Eqns. ( |38| ) for (x 2 (t)) (with 
white Gaussian noise) as a function of time in the special 
case of zero damping, i.e. 7 = 0, can be obtained by 
solving the equivalent Fokker-Planck equation [Q for the 
probability density f{x,p,t): 



^/(*,P,*) = 

d dFi (p, x) 
^ dx 



dp 



r {x) w 2 



f(x,p,t) (Al) 



The expectation value of any function M(x,p;t) can be 
written as 



(M(x,p)) 



+ 00 



dxdpM(x,p)f(x,p, t) 



(A2) 



Equations QAlJ) and ( A2) can be used to yield a BBGKY- 
like heirarchy for the evolution of phase space moments. 
Since the system we are considering is linear, this heirar- 
chy truncates exactly and yields a group of coupled lin- 
ear ordinary differential equations for the moments (x 2 ), 
{xp), and (p 2 ). These equations can be written as a single 
third-order time evolution equation for (x 2 ): 



d 3 (x 2 ) 



-477 



,d{x 2 



2a 2 (x 2 ) 



(A3) 



(A4) 



dt 3 "' dt 

subject to the initial conditions 

(x 2 (0))=x 2 (0) 
(x 2 (0)) = 2x(0)p(0) 
(x' 2 (0)) =2p 2 (0)-27]V(0) 

This equation has an analytical solution written as 

(x 2 (t)) = c x exp(r x i) + c 2 exp(r 2 i) + c 3 exp(r 3 i) (A5) 

where ci, C2, and C3 are constants depending on initial 
conditions, and n, r 2 and r 3 are the roots of a third order 
alegbraic equation 



1/3 



2a 2 - 4t7 2 x - x 3 = 

which gives 

ri = (y64/ 277/ 6 + a; 4 + a 2 ) 

- (y64/27?7 6 + a 4 -a 2 ) 
r 2 = + (V 64/277?' 

--(1 - V3i) (V64/2777 6 + a 4 + a 2 ) 



(A6) 



1/3 



,6 4_ a 4 _ a 2 



1/3 
1/3 



7*3 = r 2 



(A7) 



where the superscript * represents complex conjugation. 
The positive real root r\ implies that (x 2 (t)) will have an 
exponential growth in time. 
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